home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / MRQMIN.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  81 lines
  1. PROCEDURE mrqmin(x,y,sig: glndata; ndata: integer;
  2.        VAR a: glmma; mma: integer; lista: gllista;
  3.        mfit: integer; VAR covar,alpha: glncabynca;
  4.        nca: integer; VAR chisq,alamda: real);
  5. (* Programs using routine MRQMIN must define the types
  6. TYPE
  7.    glndata = ARRAY [1..ndata] OF real;
  8.    glmma = ARRAY [1..mma] OF real;
  9.    gllista = ARRAY [1..mma] OF integer;
  10.    glncabynca = ARRAY [1..nca,1..nca] OF real;
  11. and the variables
  12. VAR
  13.    glochisq: real;
  14.    glbeta: glmma;
  15. in the main routine. Also note that this routine calls MRQCOF, which
  16. requires a user-defined procedure FUNCS, described in that routine.   *)
  17. LABEL 99;
  18. VAR
  19.    k,kk,j,ihit: integer;
  20.    atry,da: glmma;
  21.    oneda: glncabynca;
  22. BEGIN
  23.    IF (alamda < 0.0) THEN BEGIN
  24.       kk := mfit+1;
  25.       FOR j := 1 TO mma DO BEGIN
  26.          ihit := 0;
  27.          FOR k := 1 TO mfit DO BEGIN
  28.             IF (lista[k] = j) THEN ihit := ihit+1
  29.          END;
  30.          IF (ihit = 0) THEN BEGIN
  31.             lista[kk] := j;
  32.             kk := kk+1
  33.          END ELSE IF (ihit > 1) THEN BEGIN
  34.             writeln('pause 1 in routine MRQMIN');
  35.             writeln('Improper permutation in LISTA'); readln
  36.          END
  37.       END;
  38.       IF (kk <> (mma+1)) THEN BEGIN
  39.          writeln('pause 2 in routine MRQMIN');
  40.          writeln('Improper permutation in LISTA'); readln
  41.       END;
  42.       alamda := 0.001;
  43.       mrqcof(x,y,sig,ndata,a,mma,lista,mfit,alpha,glbeta,nca,chisq);
  44.       glochisq := chisq;
  45.       FOR j := 1 TO mma DO BEGIN
  46.          atry[j] := a[j]
  47.       END
  48.    END;
  49.    FOR j := 1 TO mfit DO BEGIN
  50.       FOR k := 1 TO mfit DO BEGIN
  51.          covar[j,k] := alpha[j,k]
  52.       END;
  53.       covar[j,j] := alpha[j,j]*(1.0+alamda);
  54.       oneda[j,1] := glbeta[j]
  55.    END;
  56.    gaussj(covar,mfit,nca,oneda,1,1);
  57.    FOR j := 1 TO mfit DO da[j] := oneda[j,1];
  58.    IF (alamda = 0.0) THEN BEGIN
  59.       covsrt(covar,nca,mma,lista,mfit);
  60.       GOTO 99
  61.    END;
  62.    FOR j := 1 TO mfit DO BEGIN
  63.       atry[lista[j]] := a[lista[j]]+da[j]
  64.    END;
  65.    mrqcof(x,y,sig,ndata,atry,mma,lista,mfit,covar,da,nca,chisq);
  66.    IF (chisq < glochisq) THEN BEGIN
  67.       alamda := 0.1*alamda;
  68.       glochisq := chisq;
  69.       FOR j := 1 TO mfit DO BEGIN
  70.          FOR k := 1 TO mfit DO BEGIN
  71.             alpha[j,k] := covar[j,k]
  72.          END;
  73.          glbeta[j] := da[j];
  74.          a[lista[j]] := atry[lista[j]]
  75.       END
  76.    END ELSE BEGIN
  77.       alamda := 10.0*alamda;
  78.       chisq := glochisq
  79.    END;
  80. 99:   END;
  81.